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The navigation filter architecture successfully deployed on the Morpheus flight 
vehicle is presented. The filter was developed as a key element of the NASA Au- 
tonomous Landing and Hazard Avoidance Technology (ALHAT) project and over 
the course of 15 free fights was integrated into the Morpheus vehicle, operations, 
and flight control loop. Flight testing completed by demonstrating autonomous 
hazard detection and avoidance, integration of an altimeter, surface relative ve- 
locity (velocimeter) and hazard relative navigation (HRN) measurements into the 
onboard dual-state inertial estimator Kalman fiter software, and landing within 2 
meters of the vertical testbed GPS-based navigation solution at the safe landing 
site target. Morpheus followed a trajectory that included an ascent phase fol- 
lowed by a partial descent—to-landing, although the proposed filter architecture 
is applicable to more general planetary precision entry, descent, and landings. 
The main new contribution is the incorporation of a sophisticated hazard relative 
navigation sensor—originally intended to locate safe landing sites—into the naviga- 
tion system and employed as a navigation sensor. The formulation of a dual-state 
inertial extended Kalman filter was designed to address the precision planetary 
landing problem when viewed as a rendezvous problem with an intended landing 
site. For the required precision navigation system that is capable of navigating 
along a descent-to-landing trajectory to a precise landing, the impact of attitude 
errors on the translational state estimation are included in a fully integrated nav- 
igation structure in which translation state estimation is combined with attitude 
state estimation. The map tie errors are estimated as part of the process, thereby 
creating a dual-state filter implementation. Also, the filter is implemented using 
inertial states rather than states relative to the target. External measurements 
include altimeter, velocimeter, star camera, terrain relative navigation sensor, and 
a hazard relative navigation sensor providing information regarding hazards on a 
map generated on-the-fly. 


I. Introduction 


The extended Kalman filter has long been the workhorse of space-based navigation systems for on- 
orbit rendezvous, ranging from the early days of Apollo to the more recent era of the Space Shuttle. 
The Kalman filter is well-known and well-documented in the literature, dating from Kalman’s 
original work in 1960.! We pose the precision planetary landing problem as a rendezvous problem 
with an intended landing site. For the required precision navigation system that is capable of 
navigating along a descent-to-landing trajectory to a precise landing, the impact of attitude errors 
on the translational state estimation are included in a fully integrated navigation structure in which 
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translation state estimation is combined with attitude state estimation. In addition, the map tie 
errors are estimated as part of the process, thereby creating a dual-state filter implementation. 
Also, the filter is implemented using inertial states rather than states relative to the target. 

The problem that motivates our filter architecture approach is planetary descent-to-landing, 
however, the application of the proposed dual-state inertial extended Kalman filter is not limited 
to this class of problems. For example, this filter architecture was deployed on the Morpheus flight 
vehicle that followed a trajectory that included an ascent phase followed by a partial descent—to- 
landing. As the trajectory progresses, external measurements may include altitude (i.e, altimeter), 
relative velocity (e.g., velocimeter), attitude (e.g., star camera), position relative to a pre-loaded 
map (e.g., a terrain relative navigation (TRN) sensor), and a hazard relative navigation (HRN) sen- 
sor providing information regarding hazards to the targeting subsystem, as well as measurements 
relative to a map generated on-the-fly. An internal inertial measurement unit (IMU) provides mea- 
surements of non-gravitational accelerations and body rates. The location of the ground features 
are assumed to be corrupted by map-tie errors, all external sensors and the IMU are assumed to 
produce measurements exhibiting systematic and random errors, estimates of the initial position 
and velocity of the spacecraft are known with initial estimation errors assumed zero-mean with 
given state estimation error covariance, and uncertainty in the environment models (namely, the 
gravity) is assumed. 

Navigation algorithms based on extended Kalman filter architectures are model-based. Since 
this paper focuses on the navigation algorithms themselves, our attention is on the environment and 
sensor models employed in the EKF and not on the high-fidelity models that might be deployed in 
the various numerical simulations to model the true state history of the spacecraft and to generate 
the simulated measurements. The models of the translational and rotational dynamics of the 
lander within the EKF are key to accurately predicting the motion of the spacecraft. Furthermore, 
without an accurate measure of the non-gravitational acceleration and the spacecraft angular 
velocity as provided by an accelerometer/gyro package, the translational and rotational dynamics 
models would not be capable of supporting precision landing. The navigation algorithm must also 
possess accelerometer and gyro models. Also of interest in precision landing is the location of the 
landing site. In performing precision descent-to-landing navigation, it is necessary to be able to 
accurately predict the location of the landing site, accounting for the presence of map tie errors. In 
addition to the dynamical models describing translation and rotation of the lander, it is necessary 
to formulate the dynamics of the estimation errors. The estimation errors represent the differences 
between the true spacecraft state (position, velocity, and attitude) and the navigation provided 
estimated state. Other error dynamics associated with systematic errors are also included in the 
development. 


A. Organization of the Paper 


The organization of this paper is as follows: The descent-to-landing scenario and associated refer- 
ence frames are presented in Section II. This is followed in Section III by an overview of the basic 
structure of the extended Kalman filter (EKF). The necessary dynamics modeling, including a de- 
scription of the nonlinear equations, a treatment of the linearized error equations, and a description 
of the modeling of accelerometers, gyros, and a generic landing site are presented in Section IV. 
The proposed EKF sensor models are presented in Section VI. The models that are implemented 
in the filtering algorithms are presented in Section VII along with some notes on implementation 
issues, such as IMU thresholding and residual editing. A summary of the recommendations for 
preparation for implementation in the Morpheus lander is presented along with a description of 
the filter implementation as it was realized for the Morpheus flights is given in Section ??. 


B. Mathematical Notation 
The following notation is used throughout for variables: 


— Scalar quantities are given by non-bold lowercase and uppercase symbols, such as: a, A, 7, 
and IT. 


— Vector quantities are given by bold lowercase symbols, such as: a and +. 
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— Matrix quantities are given by bold uppercase symbols, such as: A and I. 


— Quaternions are given by bold lowercase symbols with overbars. Furthermore, the vector part 
of the quaternion is given by the same symbol in bold with no overbar and the scalar part of 
the quaternion is given by the same symbol in non-bold with no overbar. The magnitude of 
the quaternion is given the same symbol as the quaternion, but non-bold and retaining the 
overbar. 


a=| 2] and = q= |/ql| (1) 


q 


The expected value of a variable is denoted by E{-}. The vector or matrix two-norm is denoted 
by || - |]. The set of all real numbers is represented by R, such that an n-dimensional space of real 
numbers is IR” and a matrix which has n-rows and m-columns of real numbers is said to be in 
R"*™ 9 The Dirac delta is represented as 6(t — 7) and is defined via the integral 


| fse- nar = F0) (2) 


for any real-valued function f(t) that is continuous at t = 7. It is noted that the Dirac delta is 
zero for allt 4 7. The somewhat similar Kronecker delta is represented as 4;;, and is defined such 


that 
fee Tod (3) 
2) : : si 
: 0 iff 
The transpose is denoted with a superscript “T” and the inverse by a superscript “—1.” Further- 
more, the superscripts “—” and “+” represent a priori and a posteriori values (i.e. values which 


immediately precede and succeed a measurement update). A quantity which is expressed in a 
particular frame is given a lowercase symbol in the superscript to denote in which it is represented. 
The hat accent (*) is used to denote an estimated quantity, and the dot accent (° ) is used to 
denote a temporal derivative. The inner (or dot) produce is denoted by ©. The cross product 
matrix, denoted by [-x], is such that a x b = [ax] b, where a, b € R°. In terms of the components 
of a, the cross product matrix, [ax], is 


0 —a3 ag 
[ax] = a3 0 —ay : (4) 
—ag ay 0 


0 a3 «a2 ay 0 0 
fa|x]]}= |] az O a, | and [a\]J=] 0 a O |], (5) 
ag ay 0 0 0 a3 


respectively. Furthermore, note that for any a,b € R® the following properties hold 
[ax]b =—[bx]a, [a|x|]b=[b|x|Ja, and [a\]b = [b\ja. (6) 


Given the two reference frames, denoted by the a-frame and the b-frame, the transformation matrix 
which maps vectors from the a-frame to the b-frame (or equivalently expresses the relative attitude 
of the a-frame with respect to the b-frame) is given as T®. Equivalently, the quaternion which 
relates the a-frame to the b-frame is given as q?. 


II. Descent-to-Landing Scenario and Reference Frames 


A. Descent-to-Landing Scenario 


A schematic representation of the descent-to-landing scenario is shown in Fig. 1 that illustrates 
the thrust-coast-thrust engine configuration, the nominal and map tie corrupted landing site and 
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surface features, a measurement of spacecraft attitude, and a measurement of position via process- 
ing of an image of the planet surface. The descent-to-landing begins in a parking orbit about the 
planet or on a direct entry trajectory. The descent trajectory is a thrust-coast-thrust trajectory. 
The first thrust phase is a de-orbit maneuver designed to lower the orbit periapse where the landing 
sequence is initiated. Once the descent has been initiated, the spacecraft enters a coast phase until 
it reaches the orbit periapsis. After the coast phase, the braking burn occurs near periapse and 
the final descent maneuver nulls the remaining vehicle velocity to achieve a soft landing. As the 
trajectory is progresses, measurements of altitude, relative velocity, attitude (from a star camera), 
and position (from a terrain camera, which is subject to map tie error) are processed. A schematic 
representation of the descent-to-landing scenario, shown in Fig. 1, illustrates the thrust-coast- 
thrust engine configuration, the nominal and map tie corrupted landing site and surface features, 
a measurement of spacecraft attitude, and a measurement of position via processing of an image 
of the planet surface. 
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Figure 1. Schematic Descent-to-Landing Scenario 


B. Reference Frames 


Five basic classes of reference frames are considered: a planet-centered inertial frame, a planet- 
centered, planet-fixed reference frame, a planet-surface, planet-fixed reference frame, a spacecraft 
body reference frame, and various spacecraft sensor reference frames. 


1. Planet-Centered Inertial Reference Frame 


The inertial frame (denoted by 7) is taken to be the J2000 reference frame centered at the celestial 
body (at the center of mass of the body). The J2000 reference frame is defined via the FK5 star 
catalog with a standard epoch defined as 1.5 January 2000, or 12 P.M. (noon) on January 1** 2000 
in the Barycentric Dynamical Time (TDB) time scale.2 The TDB time scale is a form of atomic 
time which allows Newtonian physics, corrected for general relativity, to be followed. The J2000 
frame will serve as the base frame to which all other frames are related, which is consistent with 
the PDS Standards Reference.? The J2000 frame is known as the planet centered inertial (PCI) 
frame. 

The stellar reference frame (denoted by sr) is the frame in which the right ascension and 
declination of the set of stars is written. Typically, stellar reference frame is not aligned with 
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the J2000 reference frame, but is an inertial reference frame. The difference between the stellar 
reference frame and the J2000 inertial reference frame is a known constant transformation. 


2. Planet-Centered, Planet-Fixed Reference Frame 


There are two methods that we can utilize to compute the transformation from the planet-centered, 
planet-fixed (PCPF) reference frame (denoted by f) to the PCI reference frame, denoted by T , or 
equivalently by the quaternion ql . The first method mplements the mathematical model provided 
by the IAU/IAG/COSPAR Working Group on Cartographic Coordinates and Rotational Elements 
of the Planets and Satellites? within the EKF algorithm. The north pole of the body is defined to 
be the pole of rotation which lies on the north side of the invariable plane of the solar system and 
is specified by its right ascension and its declination. The remaining required angle specifies the 
ephemeris position of the prime meridian. The mathematical models of these angles are specified in 
Archinal, et al.? To validate the implementation, the implemented PCPF frame can be compared 
to the same frame as implemented in the SPICE™ system. For example, our analysis of the lunar- 
center, lunar-fixed (LCLF) reference frame accuracy relative to the SPICE computation showed 
that error angle between the frames is found to be on the order of XX arc sec and that unit 
directions of the inertial axes at the surface of the Moon resulted in maximum norm of the position 
differences of less than 6m. 

The second method uses an i-load initial transformation of the PCPF reference frame to the 
PCI reference frame at a given epoch (usually the start of the simulation or the initiation of 
the descent-to-landing trajectory) and then continuously updates the transformation employing a 
single rotation about the planet spin axis. This represents the simplest model for implementation 
in the EKF. In both cases, we assume that the transformation, aT (t), is known exactly—it is not a 
source of uncertainty. 


3. Planet Surface Fixed Reference Frame 


The planet surface fixed reference frame (PSF) is a reference frame (denoted by s) fixed to a point 
on the planet surface. This topocentric reference frame is selected such that the unit directions 
of the frame point in the east, north, and zenith directions. The frame is fully defined by the 
spherical latitude and longitude of the point at which the frame is attached to the surface. Since 
the latitude and longitude of the point are defined in a fixed reference system, the PSF reference 
frame is defined relative to the PCPF reference frame, denoted by T}. 


4. Spacecraft Body and Sensor Reference Frames 


The spacecraft body frame (denoted by b) and the spacecraft sensor frames are all frames which 
are fixed to and unmoving with respect to the spacecraft. The body frame is a frame which is fixed 
along axes of the structure of the spacecraft. The spacecraft body frame is depicted in Fig. 2 by 
the unit vector triad {xp, yp, Zo}. 

The spacecraft sensor frames are centered at each corresponding sensor and are allowed in- 
dependent alignments. Each sensor frame is defined by a unit vector triad, as shown in Fig. 2. 
Typically, the sensor frame +z-axis is defined to be along the bore of the sensor and the remaining 
axes span the plane perpendicular to the bore direction. Since both the spacecraft body frame and 
the spacecraft sensor frames are fixed to and unmoving with respect to the body, it is assumed 
that the relative orientation of each sensor frame with respect to the body frame is known via cal- 
ibration and testing. The sensor reference frames considered include the IMU case reference frame 
denoted by c, the velocimeter frame denoted by v, the altimeter reference frame denoted by a, the 
star camera reference frame denoted by sc, the terrain relative navigation reference frame denoted 
by t, and the hazard relative navigation reference frame denoted by h. Onboard Morpheus, the 
sensor suite did not include a star camera or a terrain relative navigation sensor, hence these are 
not illustrated in Fig. 2. A summary of the designations of the reference frames considered is given 
in Table 1. 
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Figure 2. Spacecraft Body Fixed and Sensor Fixed Reference Frames. (a) front view, (b) rear view. 


Table 1. Reference Frame Designations 


a inertial 

sr stellar 

f planet-centered, planet-fixed 
8 planet surface fixed 

b spacecraft body fixed 

Cc IMU case 

v velocimeter 

a altimeter 


SC star camera 
t terrain relative navigation 


hazard relative navigation 


III. Navigation Algorithm 


When the system dynamics and sensors can be described by linear models (possibly time- 
varying) and the measurement and process noises can be described by zero-mean, white noise 
processes, then the Kalman filter can be shown to be the optimal linear filter in terms of minimizing 
the mean-squared state estimation error.° If additionally, the measurement and process noises can 
be described by zero-mean, Gaussian white noise processes, the Kalman filter is the optimal filter. 
However, in our case, the underlying spacecraft entry dynamics and the sensors are inherently 
time-varying and nonlinear. The continuous-discrete extended Kalman filter (EKF), used as the 
primary navigation algorithm, is a direct extension of the optimal linear Kalman filter in the 
nonlinear setting (see, for example, Gelb* for more details). The EKF has long been the workhorse 
of space-based navigation systems, dating back to the early days of the Apollo missions. The EKF is 
a recursive data processing algorithm that is capable of asynchronous fusion of measurements from 
various sensors, with potentially time-varying noise characteristics. Moreover, the EKF utilizes 
prior knowledge regarding the state of the system and the error statistics associated with the state 
estimate. Additionally, the EKF is model-dependent, thus making accurate modeling of the sensors 
and the environment a key factor in the successful implementation. In order to accommodate the 
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Figure 3. Reference Frame Topology. 


nonlinearities of the system the EKF utilizes approximations based on Taylor series expansions of 
the nonlinear system dynamics and sensor models along the current estimated trajectory. With 
proper tuning and analysis, the EKF can provide a reasonably accurate estimate at each point 
along the trajectory. Tuning of the EKF can be assisted through use of Monte Carlo analysis that 
provides a means for computing the expected state estimation error covariance which can then be 
compared to the EKF state estimation error covariance. A well-tuned EKF will possess a state 
estimation error covariance that closely matches the sampled Monte Carlo covariance. 


A. EKF Architecture 


The nonlinear system model is assumed to be of the form 
x(t) = f(x(t)) + M(t)w(t), (7) 


where x(t) € R” is the state of the system, f(x(t)) € R” is the sufficiently differentiable non- 
linear system model, and w(t) € R? is a zero-mean, white-noise process with E {w(t)} = 0 and 
E{w(t)w'(r)} = Q,(t)6(t — 7) where M(t) € R"*? is the process noise mapping matrix, and 
Q,(t) = QT (t) > 0 € R®*?, Vt, is the process noise spectral density. The nonlinear measurement 
model is assumed to be of the form 

Yr = be(Xx) + (8) 


where the subscript “k” denotes a discrete time measurement at time t = tz, hy(x~) € R™ is 
the sufficiently differentiable nonlinear measurement model evaluated at the state x, = x(tx), 
and the measurement noise 7; € R™ is a zero-mean, white-noise sequence with E {7,} = 0 and 
E{n,n; } = Rzdpj where Ry = Ri > 0 c€ R™*”, Vt,. Additionally, it is assumed that the 
process noise and the measurement noise are not correlated in time, hence E {w(t)n; } = 0, Vt, tr. 
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The EKF is a continuous-discrete algorithm, as illustrated in Fig. 4. The two stages are referred 
to as propagation and update. The propagation stage continues until one or more measurements 
becomes available. Once the update is completed, the propagation begins again. The state estimate 
is denoted by x(t), the state estimation error is defined as e(t) := x(t) — x(t), and the state 
estimation error covariance is P(t) := E{e(t)e'(t)}. The state estimate and state estimation 
error covariance just prior to the update are denoted by x, , and P;,, respectively, and just after 
the update by x and P,, respectively. Note that the EKF is an unbiased filter, therefore, we 
assume that E {e(t)} = 0. For example, P; = E fe,e, °} and P} =E fefef 7}. The objective 
of the EKF is to minimize the mean-square estimation error J = Tr Py Vt. 


B. State Estimate and State Estimation Error Covariance Propagation 
Consider the propagation of state estimate and the state estimation error covariance between 
measurements during the time interval ¢ € {t,_1,t,}. The EKF propagation equations are 
x(t) = f(X(t)) (9a) 
B(t, tx-1) F(x(t)) ®(t, th) (9b) 
Q(t) = F(&))QM) + QHF* (X(4)) + M(£)Q,(t)M* (t) (9c) 


where tp_1 <t < ty, ®(t,t,_1) € R”*” is the state transition matrix mapping the state from t,_1 
to t and 


F(&(t)) := ee Paeeae Rx”, (10) 


with the initial conditions X(t,-1) = X{_,, B(te—-1,th-1) = 1, and Q(t,-1) = 0. At ty, we 
have X, = X(te), ®, = B(te,te-1), Qe = Q(t,), and P; = P(t,). After integration, the state 
estimation error covariance is mapped forward via 


P, = ©,P{_,®, + Qx. (11) 


C. State Estimate and State Estimation Error Covariance Update 
The state estimate and the state estimation error covariance update are calculated via 
Rf = &, + Ke (ye — be(X;, )) (12a) 
+ _ = aT! 
P{ = [I—K,H;,(%,)] P, [I-—K.Hie(&,)] + Ki Ri Ki (12b) 


where the Kalman gain is computed via 
Ky, = P, Hi (&,)W;": (13) 
where W,, is the measurement residual covariance matrix given by 


Wi, = Ha (& )P, Ae (&) + Re- (14) 


Altimeter, Velocimeter, Star Camera, TRN, & HRN 


Update Update Update 
—— 
IMU ——> Propagation Propagation Propagation 


Time 


I 
i ail IS, 
x, ,P 
k-1? k-1 1 gt 4 
Thy KV? Pe WP, 


Figure 4. EKF timeline. 
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and where H;(x;, ) is the measurement sensitivity matrix defined to be 


Ohy, (Xz) 


He(&) = Se e RM”, (15) 


Xk=X,, 


IV. Dynamics Modeling 


Extended Kalman filter navigation algorithms are model-based. This implies that the models 
of the translational and rotational dynamics of the lander are key to accurately predicting the 
motion of the spacecraft. Furthermore, without an accurate measure of the non-gravitational 
acceleration and the spacecraft angular velocity as provided by an accelerometer/gyro package, the 
translational and rotational dynamics models would not be capable of supporting precision landing. 
The navigation algorithm must also possess accelerometer and gyro models. Also of interest in 
precision landing is the location of the landing site. In performing precision descent-to-landing 
navigation, it is necessary to be able to accurately predict the location of the landing site, accounting 
for the presence of map tie errors. In addition to the dynamical models describing translation and 
rotation of the lander, it is necessary to formulate the dynamics of the estimation errors. The 
estimation errors represent the differences between the true spacecraft state (position, velocity, 
and attitude) and the navigation provided estimated state. Other error dynamics associated with 
systematic errors are also included in the development. 

The total acceleration of the vehicle is comprised of accelerations due to conservative forces (e.g., 
gravity), and accelerations due to non-conservative forces. The non-conservative forces acting upon 
the vehicle may include, but are not necessarily limited to, forces due to engine thrust and forces due 
to atmospheric interaction with the vehicle. An IMU consisting of a set of accelerometers and gyros 
senses the linear non-gravitational acceleration and the rotational rate of the vehicle, respectively. 
We navigate the IMU navigation base (rather than the center-of-gravity). The translational and 
rotational dynamics are 


malt) = vival) | (16a) 
mult) = Ty(tlag(rf,(t)) + TH a(d)Teas, (¢) (160) 
a) = 52h oate (160) 


where ri,,,,,(t) and vi,,,,,(£) are the position and velocity, respectively, of the IMU represented in 
the inertial reference frame, q?(t) represents the orientation of the spacecraft body fixed reference 
frame with respect to the inertial reference frame, and 


rf (t) = T/(t)r Vimu (t) a Tit (t)T 3 (G3 () 0 9 simu (t) (17a) 


fa ). (17b) 


where the non-gravitational acceleration, af,,(t), and body rate, w5,,(t), are sensed by the IMU 


9} 74 (t) 


in the IMU case reference frame, and a,(r/,(t)) is the gravitational acceleration represented in 
the planet-centered, planet-fixed reference frame, an janes (t) is the position of the center-of-gravity 
relative to the IMU represented in the body fixed reference frame. The poeetiog of the IMU relative 
to the center-of-gravity is given by a nominal position, denoted by r2, ee: ), in the body fixed 
reference frame and computed externally to the navigation system (and likely time-varying as fuel 
is expended during descent) and a deviation from the nominal position, denoted by Ar? 9 /imu (t) 
such that 


Vea = Foie ) + BE ice) (18) 
Also, we have Ti (qb(t)) = [T2(q?(t))]* where 
T? (a(t) = Isxs — 2a (t) [a2(t)x] + 2 [ab(t)x]”. (19) 


The transformation matrices Ty (t) and T¢ are assumed known. 
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A. Gravity Modeling 


There are various possibilities for modeling the gravitational acceleration, ay(rJ, (t)) in Eq. (16). 
Typically, we model the primary gravitational body using point-mass approximations and zonal 
harmonic approximations.® *:!° The simplest model assumes a spherical planet, and 


a, (rf, (t)) = - “rf 0), (20) 


where reg(t) = lez, (t)|| € R and yp is the gravitational parameter. The associated gravity gradient 
matrix is u e 
G(r, (0) = KR (3h (OL (0) = 12,()lsx3) € °°, (21) 


The gravity gradient matrix is used in propagation of the estimation error equations. 


B. Inertial Measurement Unit 
1. Accelerometers 


The IMU measurement of the non-gravitational acceleration is corrupted by errors due to 
nonorthogonality and misalignment of the axes, I, € R°®**, errors due to scale-factor uncertainties, 
S, € R®°*?, random biases, bz € R®, and noise, 7,(t) € R’. Incorporating these error sources, 
the measured non-gravitational acceleration, af. m(t) € R?, can be written in terms of the true 


non-gravitational acceleration, af,,(t) € R®, as‘)! 


Ang m(t) = (Isx3 + Ta) (Isxs + Sa) (ang(t) + ba + malt) 5 (22) 
where 
0 "V isx —~Yany Say 0 0 
Tr, = —~Yayz 0 Yaya 2 Sa = 0 Say 0 

Yazy Varn 0 0 0 Sa, 

Gm, (t) az, (t) ba. Na, (t) 

Ang mt) = | ar, (t) | »Ang(t) = | aG(t) | ba = | ba, | Malt) = | na, (t) | » 
an, (t) az (t) ba. Na. (t) 


and where 77,,(¢) is a zero-mean, white noise process with covariance 
E {na(t)na (r)} = Qad(t — 7) € R**?, 


and Q, = Q? > 0. To first-order, the non-gravitational acceleration in Eq. (26) may be written 
in terms of the measured non-gravitational acceleration and the model parameters as 


aya (t) = Ting me) _ MC awe (t))Ya = Malay om (t))Sa =i’ Ne (t) 7 (23) 
where 
| —a®,_(t) ac (t) 0 0 0 ] 
Mi (ang m(t)) = 0 0 am (t) arn, (t) 
0 0 0 0 —Am,(t) &m, (t) 

Yary 
Gn) 0 0 Yare Sa: 
Mile n@)= | 0 a) 0 |ayet= | | andsc= |, 

Ya 

0 0 at, _(t oe Sa. 

.( ) Vaz 

Yazy 
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We model ba, ¥,,, and Sz as unknown constants with 
b,=0, ¥,=0,8 =0, (24) 


and E{b,} = 0; E{y,} = 0, E{a,} =0, E{b,b]} = Pp, = PL € R**, E{a.s7} = = 
PZ, € R°%3, and E{ya7ya} = Py, =P}, € R®*®, and where Pp, > 0, Ps, > 0, and P,, > oo are 
given. In ne EKF implementation: we otiea nae the unknown constants as random constants 
driven by zero-mean, white noise processes with given process noise strength or as exponentially 
correlated random variables with given time constants and process noise strengths (see Gelb*). For 
example, the dynamics of the unknown constants in Eq. (24) can be modeled as 


ba(t) =m, (#),  Yalt) = ,, (4), Salt) = me, (0): (25) 


where E {mp, (t)} = 0, E{n,,()} = 0, E{ng, (t)} = 0, for vt, E{m, (ng, (7)} = Qv. dt — 7), 
E {n.. (t)ne, (r)} = Q,,,6(t— 7), and E {n,, (t)ny, (r)} = Q,, 6(t — 7), for Vt,7 and where Qp, = 
Qi, < R*3, Q., = QT ¢ R**?, Q,, = QT € RS, where Qn, > 0, Qs, > 0, and Q,, > 0 are 


Sa 
given. 


2. Gyros 


The IMU measurement of the angular velocity of the spacecraft is corrupted by errors due to 
nonorthogonality and misalignment of the axes, [,, € R®*3, errors due to scale-factor uncertainties, 
S, € R®*%, random biases, by € R®, and noise, 1,(¢) € R*. Incorporating these error sources, the 
measured angular velocity, wf Jism(t) € R°, can be written in terms of the true angular velocity, 


wy ,(t) € R®, as 


wEim(t) = (Isxs +p) (sxe + 8p) (wh,i(t) + by + g(t) | (26) 
where 
0 V9uz = Ving 59x 0 0 
C= 9,2 0 Vous ; So i= 0 sy 0 ; 

Yoru Vee 0 0 0 8%, 

Win, (t) ws (t) bg. Ngx (t) 

We jim(E) = | whe, ) | weet) = | w(t) | Bg = | bg, | y(t) = | M9, (t) | » 
Win, (t) w(t) bg. Ng. (t) 


and where 7,(t) is a zero-mean, white noise process with covariance 


E {ng(t)ng (7)} = Q,5(t-7) ¢ R**, 


and Q, = Qi > 0 is given. To first-order, the angular velocity in Eq. (26) may be written in terms 
of the measured angular velocity and the model parameters as 


wi i(t) = Whim (t) — Ni (W§ im (t))%q — N2(W8im(t))89 — by — 0, (t)- (27) 
where 
—tWin,(t) Wry At) 0 0 0 0 
0 0 0 we, (t) wr, () 
Voay 
et) 0 0 a Sqn 
No(w%jim(t)) = | 0 — wS,,(E) Aig | | 4-and ay 2= | ai 
0 0 wf, (t) se Sa, 
. Vex 
Vazy 
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We model bg, Y,, and sy as unknown constants with 


b,=0, 7¥,=0, and s,=0, (28) 
and E{b,} =0, E{y,} =0, E{s,} =0, and E{b,b;} = Py, = ES. eR? E{s,s)} =P,, = 
Pi €R°*S, and E We = Py, = Pe € R°*® and where Pp, > 0, Ps, > 0, and Py, > 0 are 
given. As with the accelerometer unknown constants in Eq. (24), we can model the unknown gyro 
constants in Eq. (28) as random constants driven by zero-mean, white noise processes with given 


process noise strength or as exponentially correlated random variables with given time constants 
and process noise strengths. 


C. Landing Site 
The nominal landing site position in the inertial reference frame is 
ris(t) = T} (rj, (29) 
where ri is the fixed position of the landing site in the planet-centered, planet-fixed reference frame 


and T+ (¢) is the known transformation matrix from the planet-centered, planet-fixed reference 
frame to the inertial reference frame. The planet is modeled as an oblate spheroid (or ellipsoid), 


and therefore the geoodetic coordinate system is used to represent the landing site as®§ 
(N + h) cos ¢cos » 
eee (N+h)cos¢sinA |, (30) 


(N(1—e?) +h)sing 
where h is the altitude above the reference ellipsoid, ¢ is the geoodetic latitude, is the longitude, 
e is the first eccentricity, and N is the radius of curvature in the prime vertical given by 


ye he 1) 


V1—- e2sin2 6 
where R, is the equatorial radius. While the surface features may be well-known relative to one 
another, the precise location of surface features (including the desired landing site) in the planet- 
centered, planet-fixed reference frame may possess map-tie errors, and therefore the location of the 
surface point as represented in the inertial reference frame possesses uncertainties. We model the 
true altitude, geodetic latitude, and longitude in Eq. (30) and Eq. (31) as 


h=hst+ha, ¢=¢s+¢a, and A=As +a, (32) 


where the nominal altitude above the reference ellipsoid, geoodetic latitude, and longitude of the 
landing site, denoted by his, @7s, and X75, respectively, are assumed known a priori, and the 
associated map-tie errors are denoted by ha, ¢a and Aq, respectively. We model the map-tie 
errors are unknown constants with 


ha=0, ga =0,An = 0, (33) 
where E {ha} = 0, E{¢éa} = 0, E{Aa} = 0, and E {hA} =P,, >0€ER, E{¢A} = Ps, >0€ER, 
and E{A,} = Py, > 0 € Rare given. 

D. EKF State Dynamics Model 


The state dynamics model in the navigation algorithm follows from Eq. (16)—Eq. (32). The 
EKF state dynamics are given by 


ta = 9 at) (34a) 
Vinult) = T}(t)ag (#2, (t)) + T3(G3(t)) T2AS,, (t) (34b) 
7 ile A 
q(t) = sebyilt) @ Gi (t) (34c) 
ba = 0,44 = 0,8 = 0, by = 0,4, = 0,8) =0,ha =0,¢a =0,An =0, (34d) 
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where 


eG) = TiOPmulé) + TE OTH), jmu(t) (35a) 
TG) = Isx3 — 242(¢) [ab(e)x] + 2 [42() x]? (35c) 
“ Teas i (t 
Wp i(t) = ( ( ' (35d) 
and 
ang (t) = oe (t) - M, (Ay aaal lq — M2 (@),9,m)8a 7 hy (36a) 
Eel) = WE gn 8) — Ni (8 an) ¥q — Na(W8 ign 8g — Bg (360) 


V. Estimation Error Dynamics 


A. Attitude Estimation Errors 
Define the attitude estimation error as 
= _ a —1 
day (t) = a(t) @ [a] (37) 


Taking the time derivative of the error quaternion in Eq. (37), substituting q?(t) and q(t) from 


Eq. (16c) and Eq. (34c) respectively, using the definition of quaternion multiplication, and re- 
arranging yields 


ST as 


—dw},(t) © daz (t) 
where dw} ),(t) == wp,,(t) — cp; (t) =—T" (w5), (t) — w(t). Assuming small attitude errors and 
neglecting second-order terms, Eq. (38) reduces to 

: 7 1 
day(t) = —@%,,(t) x dah(t) + 5 oil) (39) 
dqz(t) = 0. 


Define dby := by — By OSq ‘= Sg — 8g, and dy, := Y, — Yg. Computing dw ;(t) using Eq. (27) 
and Eq. (36b) yields 


Sup; (t) = -T? [Nu W5 im (O)O%g + N2(W5/im(t))68q + dbg + ng(t)| . (40) 
Substituting Eq. (40) into Eq. (39) yields 
é Cc 1 Cc 
Bab (t) = = [Tw Es m())] dab) — STE (Na Fyim (OEY (41) 


+ Na (10 jim (€))58y + Sby + 14(t)) - 
Assuming small attitude error angles, denoted by da, we have the approximation 
dq? (t) & 55a). (42) 
Taking the time-derivative of Eq. (42) and substituting Eq. (41) yields 
B6x(t) = — | Trer§ 5m] See(t) — TE (Na (5 jen ())O%y + No Wh pim(t))I89 + dy + g(t) - (43) 


With the attitude estimation error as in Eq. (37) and the approximation in Eq. (42), we have the 
relationship 


Ti(at(O) = TYG) | + (Sar |. (44) 
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B. Position and Velocity Estimation Errors 


Define the position estimation error and velocity estimation as 6rt,,,,,(t) = Lm (t) — Pima (t) 


and 6viny(t) = Ving (t) — Vin, (t), respectively. Taking the time-derivative, using the relationships 
in Eq. (16) and Eq. (34), and collecting terms yields 


SFimu(t) = SVimu(t) (45a) 
SVimu(t) = T(t) (ag(rég(t)) — ag (#Z,(t))) + Th (G2 (t)) Tear g(t) — Ty (Gi (t))Teaig(t). (45b) 


Expanding the acceleration due to gravity in a Taylor series and keeping the first-order terms 
yields the approximation 


a,(r2q(t)) — ag (FZ, (t)) + G(#,(t))oré,(t), (46) 


where the gravity gradient matrix G(r/,(t)) is given in Eq. (21) and where we define dr{,(t) := 
r/,(t) — #4,(t). Note that we implicitly assume that the true gravity has the same functional form 
as the model employed in the EKF-in this case, a spherical gravity model. In reality, the true 
gravity is significantly more complex; however, this approach has proven to be effective in our 
EKF filter implementation. More complex gravity models can certainly be employed with the 
associated additional complexity in the state and state estimation error covariance propagation. 


Using Eq. (17), Eq. (18), Eq. (35a), and the definition OD Hien) = Art jimut) - AF a imult), 
it can be shown that 
rl, (t) = Ti ()orinu(t) — TE OTBO) [Pea pimu] 5a) + TAO) jimall) (47) 
Then from Eq. (27), Eq. (36), Eq. (44) it follows that to first-order we have 
Ty (G2 (t)) Tear g(t) — TG?) Teang(t) = —Ty (4?) [Taig m(t) x] Ser(t) (48) 


~T (G7 (t)) Te (Sba(t) — Mi (@hig,m (t))5-¥a(t) — Maing m (t))58a(t) — a(t) 
— GEL, ()) TY (GF) ore 9 imu (t), 


where dbg = ba — ba, 08a = Sa — 8a, and OY, = Yq — Yq. Substituting Eq. (47) and Eq. (48) into 
Eq. (45b) yields 


— [T4(43O) [Teak m(t)>] + GL (O)TH EO) [Fy simu] ] Fert) (49) 
—T),(G2(1))T? (Sb a(t) + Ma (Aig, (t)) Ya (l) + Ma (Align (€))58a(t) + M4 (0)): 


VI. Sensor Modeling 


The EKF is a model-based algorithm requiring models of the sensors. In general, there are two 
classes of sensor models required for high-fidelity analysis of the integrated guidance, navigation, 
and control (IGN&C) system. These are (1) high-fidelity models to support the high-fidelity 
analysis in simulation, and (2) sensor models for the navigation algorithm. This paper is concerned 
only with sensor models for the navigation algorithm. To that end, a navigation sensor model 
includes the following: 


1. A mathematical model represented by a nonlinear equation(s) as a function of the states of 
the system (such as position, velocity, attitude, map tie, etc.). 


2. A measurement mapping matrix comprised of the partial derivatives (of the model above) 
evaluated at the most recent state estimate. 


3. An error model comprised of random noise and systematic errors, including representative 
values of the uncertainty in the various error sources. 


The sensors currently under consideration are: 


14 of 24 


American Institute of Aeronautics and Astronautics 


1. Altimeter 
2. Velocimeter 
Star camera 


Terrain relative navigation (TRN) 


Cr, aie es 


Hazard relative navigation (HRN) 


The models presented here represent the class of sensors that might be used to support precision 
EDL. Not all sensors will likely be used in all phases of EDL. For example, for the Morpheus 
project, the star camera and terrain camera (TRN) were not employed at any time. Additionally, 
the models used in the EKF are dependent on the actual sensors selected. For example, there are 
various possibilities for the landing altimeter (first return, slant range, etc.). This paper presents 
an illustrative set of models. 


A. Altimeter 


The altimeter is modeled as 
Nsphik = tex — Vgrd,k a batt - Talt,k (50) 


a es i(gb )pb 
where Tar. ¢ = |lPimu,k + T5(Gi,¢)Pait/imulls Pora,k is the distance of the sub-spacecraft surface point 


from the center of the planet, r is assumed known, bai is the altimeter sensor bias, and 


b 
alt/imu 
Nalt,k iS a Zero-mean, white noise sequence with E {nai,x} = 0, Vt, and E {nine} = Ratt,nOn,j € R, 
Rait,n > 0, Vtz. We often model the altimeter bias, bai, as a random constant with 


bat = 0, (51) 


where E {bai} = 0 and E VOcer = P,,,, > 0 €R. As with other biases (e.g., IMU), we might employ 
a random constants model driven by zero-mean, white noise processes with given process noise 
strength or as exponentially correlated random variables with given time constants and process 
noise strengths (see Gelb*). The are various ways to model rgra,~ in the EKF. For example, we 
might model the planet as a sphere and rgrq,z would be the spherical radius of the planet plus a 
term to account for the local surface topography (possibly the mean height of the planet surface 
above the spherical radius) plus an additional time-varying term to account for deviations from 
the mean height and for small to medium surface features (e.g., rocks and small craters). 
The associated estimate of the altimeter measurement is computed via 


Asph,k = Palt,k P ord,k “++ bait (52) 


si Jaa i(Gb pb * 
where Pip ¢ = llPimue + Tb (G2,4)¥ait/imull> Pgrd,k is an estimate of the distance of the sub-spacecraft 


surface point from the center of the planet, and Deis is the altimeter sensor bias estimate. Defining 


ONsph,k = Nspn,k — Nspn,k, it follows that to first-order we have 


Palt,k we 
ONsph,k — ( re a Ts (42g) [ote simu x]Baxx) a Or grd,k + Obait + Nalt,k > (53) 
alt,k 


qt 


Aq ay ifadb b ce r oe A 
where Tale — Tick + i (G5) Fait /ienu? Obait — batt batt and Or grd,k “= Vord,k — Vgrd,k- 


B. Velocimeter 


The velocimeter is modeled via 


Vee i = T?T3 (a? ,) [Yet wri k x tae + Dyer + Nyelk (54) 


15 of 24 


American Institute of Aeronautics and Astronautics 


where 
Yep =r + Tig? .\e° ves. P(g’ (wn ee (55) 
vel,k imu,k b Dik vel/imu? ‘vel,k imu,k b ik b/i,k vel/imu } ? 


b 
Vvel/imu 
fixed reference frame of the planet with respect to the inertial reference frame represented in the 


inertial reference frame and wi, i is the constant known rotation vector of the planet about the spin 


is assumed known, wy c= Ty. nw i is the angular velocity of the planet-centered, planet- 


axis, r’., , and Me , are the position and velocity of the velocimeter in the inertial reference frame, 
respectively, bye; is the velocimeter bias in the velocimeter reference frame, and 1,,,; ;, is a zero- 


mean, Gaussian, white-noise sequence with E {Nvet,k } = 0, Vt, and E { rivet a Mere} = Roet,k or, € 
R°*3, Ruei,~ > 0, Vtx. We often model the velocimeter bias, byi, a8 a random constant with 


bye = 0, (56) 


where E {buyer} = 0 and E{ byeib2.,} = Pp 


measurement is 


> 0 € R®*%. The estimated relative velocity 


vel 


aU __ mumi/sb a4 a at ‘f 
Vrel,k ~_ T;T3(G: x) (Wert _ Wik x fae alr Duel, (57) 
where 
at _ ai 4 (2b b a4 — at i (2b ab b 
Poel,k ~ Vimu,k a TAG: gt ier Vvel,k ~~ Vimu,k a Ti (Gi.x) (Bhi x i pe) ’ (58) 


With OV rel, k = Verel,k — Vrel,kr ODvet = Buet — Byer, and Eq. (40), it follows that to first-order we 
have 


OViel,k =~ T?T; (a? ,) (Ges OF ia _ SVE nu) a :| [ (2b isms x tt ipeeh) x] 


? 


+ TG) [oot nr] THA) [roesyimuX] + [THA e) (Phar — Oye % Ber) *] |>0 (59) 
+ TT; pee x] T’ (ob, + N1(W5 ji mg) Vg + N2(5 sim, 589 + na) + dbvet + Net, 
where Wh ime is the measured angular velocity at t, and WE fim = Tee img" 
C. Star Camera 
The star camera is modeled via 
Gek = Iy.k @ 5° @ Gi, @ Gr, (60) 


where q#° and qi,. are known transformations, and q,,; is the bias-noise quaternion given by 


i,k = > 


where 0% = Bsc + Nsc,k1 Ok = ||@x||, and n,.., is a zero-mean, Gaussian, white-noise sequence with 
E{n,.,} = 0, Vty and B {necaMne} = Rocpdkj € R33, Rec. > 0, Vez. We often model the 


star camera bias, b,-, as a random constant with 


bs. = 0, (61) 


where E {b,-} = 0 and E {b,-by.} = Pp,, > 0 € R®*°. 
The estimated star camera quaternion is 


Oo . = Gna GS al, @ ai. (62) 
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where 0; = bsc. With the definitions 5q%¢, = 4%, ®@ 450, , 642, = 4, @G, , and 5Gy,4 = 
Gr,k ® Te and from Eq. (60) and Eq. (62), it follows that to first-order we have 


5n,& + T(Gn,«)T§5a? , 


6Gsrk = (63) 
1 
Assuming small angles and with the definition dbs. := Dse — Dies we have 
1 1 
OQn,k = 3 oP se - 3 Isc,k (64) 
Then with dq?, + 50a, and dq%* ,, © 56a, substituting Eq. (64) into Eq. (63) yields 
ow = T (Gyn) TZ oar + dbse + Neck (65) 


D. Hazard Relative Navigation 


HRN tracking is comprised of an hazard detection and avoidance (HDA) phase to create a 
feature map followed by a navigation phase in which HRN measurements are processed. The 
feature map specifies the locations of the key features relative to the corrected intended landing 
site (ILS) where the corrected ILS incorporates that latest estimates of the map-tie errors. As 
the descent progresses and when the resolution dictates, a new feature map is generated—this 
mapping process is invisible to the navigation algorithm except that HRN measurements are not 
provided for short periods during the mapping—and the navigation phase continues once the HRN 
measurements are again provided. During each navigation phase, we assume that only one feature 
is tracked. 

The first step in the HDA process is for navigation to provide to the HDA algorithm the current 
inertial state estimate for the purpose of locating the HRN sensor with respect to the corrected 
ILS. The location of the HRN sensor relative to the corrected ILS is given in the planet-surface, 
planet-fixed reference frame and the location of the feature relative to the HRN sensor is given in 
the HRN sensor reference frame. 

The corrected ILS, represented by Te in Fig. 6, is given by Eq. (30) where his, dis, and Ais, 
respectively, are assumed known a priori, and the associated map-tie errors denoted by ha, ba 
and Ads respectively, are estimated before the HDA process starts (likely during the TRN phase). 
As depicted in Fig. 6, the estimated location of the HRN sensor relative to the corrected ILS at 
t=t, is given by 


PTLS/HRN,j = —T¥ eee (Fars + Ti )rew/rmu) = itis) , (66) 
where 
—sin(Ais + Aa) cos(Ais + Aa) 0 
* = | —sin(¢is + Ga) cos(Ars + Aa) sin(dis + da) sin(Ais + AA) cos(dis + Pa) | ; 


cos(¢is + Ga) cos(Mis +Aa) —cos(dis + Ga) sin(Ais + Aa) — sin(dis + da) 


and ty RN/IMU is known. During the mapping phase, we compute the locations of tracking features 
in the planet-surface, planet-fixed reference frame with the origin at the corrected ILS. Referring 
to Fig. 7, we obtain the relationship Pr/ILs as 


« 7 i (2d beh Rn 
Pr/ILS = 8 Co gm re _ PTLS/HRN,m ’ (67) 


where t = tm, ? an and fing /HRN,m 2te propagated from t = t; to t = tm using an HRN-based 
IMU initialized by state estimates from the EKF to compute a ; and firs /HRN,j? respectively, 
and v5 /HRN,m is the HRN measurement where, 


Yb/HRN.m = Vb/HRNm + Darn + Qhenym (68) 
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Truth Navigation 
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Soe a nk SS 7 See 


Figure 5. HDA modeling—Step 1—Locating the HRN sensor. 


Darn is the HRN bias in the HRN reference frame, and 7),.,, ,, is a zero-mean, Gaussian, white-noise 


sequence with E {tian = 0, Vim and E { thienm"ben a} ~ Rarnymdmj € Ree, Rairn,m > 0, 
Vtm. We often model the HRN bias, by;n, as a random constant with 


Darn = 0, (69) 


where E {barn} = 0 and E {Darnbien } = Phim > 0 € R°*3. 


Arn 


Consider the navigation phase of the HRN depicted in Fig. 8. Using the measurement model 
in Eq. (76) we can write 


oh oh w 
Yr/HRNym = 'F/HRN,m + Phrn- (70) 
To compute by by RNm °= Yb/HRN.m ~ VF/HRN m: Consider Fig. 8 to obtain 
h _ ns TT qo a b _T qb TT af pf aS —A s 
lr/HRN,k b | Ls(Gie)Prmu,e + YHRN/IMU b (Gin) pk \Trns + 1s \Tryits lpr , 
from which it follows that 


PP / HRN. =-T} Tila dtc. + riERN/IMU _ Ts (Gi x) T 7 (th) (tts + Tr (fais = Ari.x))| : 


Defining dbprn,k = bina ~biee® JATE = Ari, —Ar%,;, and Or ey HRN. = Th /HRN & ~ EP) HRN, k? 
it follows that 
SY es HRN m = Orb / HRNK + 6Darn,k + Mhrn,k> (71) 
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Truth Navigation 
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=a) eo Se 
Figure 6. HDA modeling—Step 2—Locating the surface features. 
where 
h Ari (Ab a Ari (2b io s 
Ores HRN, k =-|T; Ti (G:.x) Ort MU, k = 5, T3(a? ,)T+,T! dAL® (72) 


=i Ti (G21) (as ~ Thr (fas +05 (Fyre ~ Ari.) ) x |dQp. 


VII. Navigation Algorithm Implementation 


Let the state vector, x, at t, be given by 
So cg/imu | ¢ p49x49 


(73) 
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Figure 7. HRN navigation phase. 


where the parameter vectors for the accelerometer, gyro, map tie errors, and sensor biases are 
defined as 


ba 
ba by ha - . 
bo |. gy |e DeH |S, la Dae=| oa |y o0d pe= a 
r sc 
Ya V9 . birn 


Therefore, pa € R!’, p, € R'”, pm € R®, and p, € R'™. 

The TRN and HRN will not be used for navigation at the same time. This enables the reuse 
of several key state variables. During the TRN phase we expect that the map-tie errors will 
be reduced, hence we will estimate pm = (ha, da, Aa). Once the HRN phase begins, we stop 
updating the corrected ILS to ensure that guidance and targeting are not chasing a moving landing 
site. However, during HRN we will need to estimate Ar*, € R®, so we reset p,, > Ar%. We can 
also reset Dipn 4 Darn. 

One caveat to the state estimate update given by Eq. (12) is that of the attitude update. All of 
the preceding developments have focused on the attitude error in the form of small angles, denoted 
by da,. If the portion of 7 pertaining to the attitude is given by 6a; , then the quaternion 
update is 
® q! kp (74) 


Gik = 


Lea t 
Abt 5 0G, 


where q?;, is the a priori estimate of the quaternion and 6@, = 0. The quaternion in Eq. (82) 
satisfies the unity norm constraint only to first order, and so to ensure that the quaternion remains 
unity norm a normalization is performed. 
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